randn('state',32);
rand('state',32);
load oillogitdata;

id=data(:,1);
s_jt=data(:,2);

x1=[data(:,3:33)];
clear data;
load oilblpiv;
iv=data(:,2:24);
brand=floor(id/100000);
company=ones(1680,1);
pre_merger_dum=cr_dum(company);
own_dum=pre_merger_dum(1:7,:);



load oiloutside;
meanprice=data(:,1);
meanshare=data(:,2);
meany=log(reshape(meanshare,7,10))-log(ones(7,10)-kron(ones(7,1),sum(reshape(meanshare,7,10))));


%IV has columns 1 to 9 containing the average price across each region excluding that corresponding to the current observation in each of the 9 time periods.  
%Columns 10 to 15 contain brand dummies
%Column 16 to 24 contain time dummies

%IV=[iv(:,2:2) x1(:,2:14)];
IV=[iv x1(:,4:31)];
clear iv data;
horz=['OLS    IV'];
vert=['constant  ';
      'price     ';
      'private   ';];
%10 REGIONS...7 brands...24 months...
nmkt =24*10;     % number of markets = (# of cities)*(# of months)  %
nbrn =7;     % number of brands per market. if the numebr differs by market this requires some "accounting" vector %
cdid=kron([1:nmkt]',ones(nbrn,1));
cdindex = [nbrn:nbrn:nbrn*nmkt]';       
% create weight matrix
invA = inv([IV'*IV]);



% compute the outside good market share by market
temp = cumsum(s_jt);
sum1 = temp(cdindex,:);
sum1(2:size(sum1,1),:) = diff(sum1);
outshr = 1.0 - sum1(cdid,:);
outshr=outshr.*(outshr>=0);
y = log(s_jt) - log(outshr);

%grand means for elasticity calculations
mean_share_atavg=mean(reshape(s_jt,7,24*10)')';
mean_price_atavg=mean(reshape(x1(:,1),7,24*10)')';

%ols logit
beta_ols = inv(x1'*x1)*x1'*y;
res = y - x1*beta_ols;
%vcov_ols = 1/size(x1,1)*res'*res*inv(x1'*x1);
vcov_ols = inv(x1'*x1)*(x1'.*(ones(size(x1,2),1)*res')*(res*ones(size(x1,2),1)'.*x1))*inv(x1'*x1);
se_vcov_ols=sqrt(diag(vcov_ols));

%Newey S.E.'s
%Need to get lagged value of xe

x1_lag1=x1(1:1680-70,:);
x1_lag2=x1(1:1680-70*2,:);
x1_lag3=x1(1:1680-70*3,:);


x1_temp1=x1(70+1:1680,:);
x1_temp2=x1(70*2+1:1680,:);
x1_temp3=x1(70*3+1:1680,:);


res_lag1=res(1:1680-70,:);
res_lag2=res(1:1680-70*2,:);
res_lag3=res(1:1680-70*3,:);

res_temp1=res(70+1:1680,:);
res_temp2=res(70*2+1:1680,:);
res_temp3=res(70*3+1:1680,:);


G_0=x1'.*(ones(size(x1,2),1)*res')*(res*ones(size(x1,2),1)'.*x1);
G_1=x1_lag1'.*(ones(size(x1_lag1,2),1)*res_lag1')*(res_temp1*ones(size(x1_temp1,2),1)'.*x1_temp1);
%Need to get twice lagged value of xe
G_2=x1_lag2'.*(ones(size(x1_lag2,2),1)*res_lag2')*(res_temp2*ones(size(x1_temp2,2),1)'.*x1_temp2);
%Need to get three times lagged value of xe
G_3=x1_lag3'.*(ones(size(x1_lag3,2),1)*res_lag3')*(res_temp3*ones(size(x1_temp3,2),1)'.*x1_temp3);



vcov_ols_neweywest = inv(x1'*x1)*(G_0+.5*(G_1+G_1'))*inv(x1'*x1);
se_vcov_ols_newey=sqrt(diag(vcov_ols_neweywest));
clear x1_temp1 x1_temp2 x1_temp3 res_lag1 res_lag2 res_lag3 res_temp1 res_temp2 res_temp3 G_0 G_1 G_2 G_3;



%ols elasticities at grand means
alpha=beta_ols(1);
o1 = alpha*mean_share_atavg*mean_share_atavg';
o2 = diag(alpha*mean_share_atavg);
omega=(o2-o1);
olselas=omega.*(1./mean_share_atavg*mean_price_atavg');




% IV regressions
beta_IV = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*IV'*y;
res = y - x1*beta_IV;
IVres = IV.*(res*ones(1,size(IV,2)));
b = IVres'*IVres;
vcov_IV = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*b*invA*IV'*x1*inv(x1'*IV*invA*IV'*x1);
se_vcov_IV=sqrt(diag(vcov_IV));


%Newey West S.E.'s here

IV_lag1=IV(1:1680-70,:);
IV_lag2=IV(1:1680-70*2,:);
IV_lag3=IV(1:1680-70*3,:);


IV_temp1=IV(70+1:1680,:);
IV_temp2=IV(70+1:1680,:);
IV_temp3=IV(70+1:1680,:);


res_lag1=res(1:1680-70,:);
res_lag2=res(1:1680-70*2,:);
res_lag3=res(1:1680-70*3,:);

res_temp1=res(70+1:1680,:);
res_temp2=res(70+1:1680,:);
res_temp3=res(70+1:1680,:);

IVres_lag1 = IV.*(res*ones(1,size(IV,2)));
IVres_lag2 = IV.*(res*ones(1,size(IV,2)));
IVres_lag3=IV_lag3.*(res_lag3*ones(1,size(IV_lag3,2)));

G_0=IVres'*IVres;
G_1=IV_lag1'.*(ones(size(IV_lag1,2),1)*res_lag1')*(res_temp1*ones(size(IV_temp1,2),1)'.*IV_temp1);
vcov_IV_neweywest = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*(G_0+.5*(G_1+G_1'))*invA*IV'*x1*inv(x1'*IV*invA*IV'*x1);
se_vcov_IV_neweywest=sqrt(diag(vcov_IV_neweywest));

clear IV_lag1 IV_lag2 IV_lag3 IV_temp1 IV_temp2 IV_temp3 res_lag1 res_lag2 res_lag3 res_temp1 res_temp2 res_temp3 IVres_lag1 IVres_lag2 IVres_lag3;



%iv elasticities at grand means
alpha=beta_IV(1);
o1 = alpha*mean_share_atavg*mean_share_atavg';
o2 = diag(alpha*mean_share_atavg);
omega=(o2-o1);
ivelas=omega.*(1./mean_share_atavg*mean_price_atavg');



%Approximate first-stage F-stat, leave out vcov d.f. adjustment as N is large relative to number of regresors...
% b_1st=inv(IV'*IV)*IV'*x1(:,1);
% res=x1(:,1)-IV*b_1st;
% vcov_1st=inv(IV'*IV)*(IV'.*(ones(size(IV,2),1)*res')*(res*ones(size(IV,2),1)'.*IV))*inv(IV'*IV);
% se_vcov_1st=sqrt(diag(vcov_1st));
% R=[eye(9) zeros(9,13)];
% q=9;
% F=(R*b_1st)'*inv(R*vcov_1st*R')*(R*b_1st)/q;

%Test
%This evaluates first order conditions at post merger averages.  This is
%done by setting the tolerance on fsolve to be extremely high and then
%looking at the function evaluation at the "optimum".  

load 'c:\data\restat\analysis_data\postoiloutside';
post_mean_price=data(:,1);
post_mean_share=data(:,2);
post_mean_y=log(reshape(post_mean_share,7,10))-log(ones(7,10)-kron(ones(7,1),sum(reshape(post_mean_share,7,10))));
IV_evaluated=zeros(10,7);
alpha=beta_IV(1);
pre_merger_dum=cr_dum(company);

for t=1:10
        shares=meanshare((t-1)*7+1:t*7);
 	    mp=meanprice((t-1)*7+1:t*7);
        own_dum=pre_merger_dum((t-1)*7+1:t*7,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*7+1:t*7,:)=o2-o1;
        elas((t-1)*7+1:t*7,:)=omega((t-1)*7+1:t*7,:).*(1./shares*mp');
        markup((t-1)*7+1:t*7,1)=-inv(omega((t-1)*7+1:t*7,:).*(own_dum*own_dum'))*shares;
end

mc_hat=meanprice-markup;
company_post=company.*(1-(company==4))+(company==4)*6;
post_dum=cr_dum(company_post);

first=1;
for i=1:10
	last=i*7;
    expall=exp(post_mean_y(:,i)-alpha*post_mean_price(first:last));
	IV_evaluated(i,:)=feval(@logitbert_eq,post_mean_price(first:last),expall, alpha, mc_hat(first:last), post_dum(first:last,:))';   
	first=last+1;
end
IV_FOCs=mean(IV_evaluated((sum(reshape(meanshare,7,10))<1)',:));


ols_evaluated=zeros(10,7);
alpha=beta_ols(1);
pre_merger_dum=cr_dum(company);

for t=1:10
        shares=meanshare((t-1)*7+1:t*7);
 	    mp=meanprice((t-1)*7+1:t*7);
        own_dum=pre_merger_dum((t-1)*7+1:t*7,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*7+1:t*7,:)=o2-o1;
        elas((t-1)*7+1:t*7,:)=omega((t-1)*7+1:t*7,:).*(1./shares*mp');
        markup((t-1)*7+1:t*7,1)=-inv(omega((t-1)*7+1:t*7,:).*(own_dum*own_dum'))*shares;
end

mc_hat=meanprice-markup;
company_post=company.*(1-(company==4))+(company==4)*6;
post_dum=cr_dum(company_post);

first=1;
for i=1:10
	last=i*7;
    expall=exp(post_mean_y(:,i)-alpha*post_mean_price(first:last));
    ols_evaluated(i,:)=feval(@logitbert_eq,post_mean_price(first:last),expall, alpha, mc_hat(first:last), post_dum(first:last,:))';
    first=last+1;
end
OLS_FOCs=mean(ols_evaluated((sum(reshape(meanshare,7,10))<1)',:))



%number of bootstrap draws
sims=1000;
bs_ols_evaluated=zeros(sims,7);

for i=1:sims
    bs_beta_ols=beta_ols+chol(vcov_ols_neweywest)'*randn(size(beta_ols,1),1);
    alpha=bs_beta_ols(1);

    for t=1:10
        shares=meanshare((t-1)*7+1:t*7);
     	mp=meanprice((t-1)*7+1:t*7);
        own_dum=pre_merger_dum((t-1)*7+1:t*7,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*7+1:t*7,:)=o2-o1;
        markup((t-1)*7+1:t*7,1)=-inv(omega((t-1)*7+1:t*7,:).*(own_dum*own_dum'))*shares;
    end

    mc_hat=meanprice-markup;
    company_post=company.*(1-(company==4))+(company==4)*6;
    post_dum=cr_dum(company_post);

    first=1;
    for j=1:10
        last=j*7;
        expall=exp(post_mean_y(:,j)-alpha*post_mean_price(first:last));
        ols_evaluated(j,:)=feval(@logitbert_eq,post_mean_price(first:last),expall, alpha, mc_hat(first:last), post_dum(first:last,:))';
        first=last+1;
    end
    bs_ols_evaluated(i,:)=mean(ols_evaluated((sum(reshape(meanshare,7,10))<1)',:));   
end

bs_IV_evaluated=zeros(sims,7);

for i=1:sims
    bs_beta_IV=beta_IV+chol(vcov_IV_neweywest)'*randn(size(beta_IV,1),1);
    alpha=bs_beta_IV(1);

    for t=1:10
        shares=meanshare((t-1)*7+1:t*7);
     	mp=meanprice((t-1)*7+1:t*7);
        own_dum=pre_merger_dum((t-1)*7+1:t*7,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*7+1:t*7,:)=o2-o1;
        markup((t-1)*7+1:t*7,1)=-inv(omega((t-1)*7+1:t*7,:).*(own_dum*own_dum'))*shares;
    end

    mc_hat=meanprice-markup;
    company_post=company.*(1-(company==4))+(company==4)*6;
    post_dum=cr_dum(company_post);

    first=1;
    for j=1:10
        last=j*7;
        expall=exp(post_mean_y(:,j)-alpha*post_mean_price(first:last));
        IV_evaluated(j,:)=feval(@logitbert_eq,post_mean_price(first:last),expall, alpha, mc_hat(first:last), post_dum(first:last,:))';
        first=last+1;
    end
    bs_IV_evaluated(i,:)=mean(IV_evaluated((sum(reshape(meanshare,7,10))<1)',:));   
end


diary output_logitfoctest_oil_coll
    'ols'
    OLS_FOCs
    'ols confidence interval'
    quantile(bs_ols_evaluated,.05)
    quantile(bs_ols_evaluated,.95)
    'IV'
    IV_FOCs
    'IV confidence interval'
    1000*quantile(bs_IV_evaluated,.05)
    1000*quantile(bs_IV_evaluated,.95)
diary off
